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ABSTRACT 

We present results on spectral index distributions, number counts, redshift distribution 
and oth er general statistical proper ties of extragalactic point sources in the NEWPSstj 
sample ([Lopez- Caniego et a,l.ll200~7h . The flux calibrations at all the WMAP channels 
have been reassessed both by comparison with ground based observations and through 
estimates of the effective beam areas. The two methods yield consistent statistical 
correction factors. A search of the NASA Extragalactic Database (NED) has yielded 
optical identifications for ~ 89% of sources in the complete sub-sample of 252 sources 
with S/N ^ 5 and S ^ 1.1 Jy at 23 GHz; 5 sources turned out to be Galactic and 
were removed. The NED also yielded redshifts for ~ 92% of the extragalactic sources 
at |&| > 10°. Their distribution was compared with model predictions; the agreement 
is generally good but a possible discrepancy is noted. Using the 5 GHz fluxes from 
the GB6 or PMN surveys, we find that ~ 76% of the 191 extragalactic sources with 
S23GHZ > 1-3 Jy can be classified as flat-spectrum sources between 5 and 23 GHz. 
A spectral steepening is observed at higher frequencies: only 59% of our sources are 
still flat-spectrum sources between 23 and 61 GHz and the average spectral indexes 
steepen from (af) = 0.01 ± 0.03 to {a\\) = 0.37 ± 0.03. We think, however, that the 
difference may be due to a selection effect. The source number counts have a close to 
Euclidean slope an d are in good a g reeme nt with the predictions of the cosmological 
evolution model by iDe Zotti et all (2005). The observed spectral index distributions 
were exploited to get model-independent extrapolations of counts to higher frequen- 
cies. The risks of such operations are discussed and reasons of discrepancies with other 
recent estimates are clarified. 

Key words: surveys - galaxies: active - cosmic microwave background - radio con- 
tinuum: galaxies - radio continuum: general. 



1 INTRODUCTION 

The statistical properties of extragalactic sources above 
~ 10 GHz are still largely unknown. This is due to the 
great difficulty of carrying out extensive surveys at high 
radio frequencies, because of the very small fields of view 
of ground-based radio telescopes and of the need of rel- 
atively long on-source integrations. Up to the end of the 
last century, the highest radio frequency at which large- 
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area sky surv e ys had been p erformed was 5 GHz (GB6 , 
Gregory et all dl996h: P MN. iGriffith et all l|l994 Il995l b 
Wright et all l|l994ll996h ). Recently, the situation has con- 



siderably improved thanks to the Pilot ATCA survey that 
covered the declination regi on between —60° and —70° at 
18 G Hz down to 100 mjy (|Ricci et all 12004 ISadler et ail 
l2006t) , to the 9C survey covering 520 deg 2 at 15.2 GHz to 
a flux limit of 25 nijy but going down to 10 mjy on small 
areas (|Waldram et al.ll2003l ; iTavlor et al.ll200ll ). and to data 
from the VSA in the dec lination range —5° < 5 < 60° at 33 
GHz (|Clearv et alj|2005h . The on-going Australia Telescope 
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WMAP band & 
central freq. (GHz) 



Average 
effective area (sr) 



Symmetrized area 
(sr) 



Self-calibration correction 
average value + error 



K (23) 
Ka (33) 
Q (41) 
V (61) 
W (94) 



(2.595 ±0.014) x 10~ 4 
(1.569 ±0.012) x 10" 4 
(1.016 ±0.007) x 10~ 4 
(4.81 ±0.07) x 10~ 5 
(3.36 ±0.17) x 10~ 5 



2.46 x 10~ 4 
1.44 x 10" 4 
8.94 x 10~ 5 
4.19 x 10~ 5 
2.07 x 10~ 5 



1.050 ±0.006 
1.086 ± 0.008 
1.136 ±0.008 
1.15 ±0.02 
(1.62 ±0.08) 



Calibrators correction 
average value + error 



0.99 ± 0.06 
1.12 ± 0.07 
1.15 ± 0.08 
1.32 ± 0.09 
(1.41 ± 0.10) 



Table 1. Results of the self-calibration process (effective area/symmetrized area; see § 12.31 1 and of the source calibration by comparison 
with ground based measurements (ground based fluxes/NEWPS fluxes or average ratio of the corrected flux to the uncorrected one; see 
£j 12.21 1. The correction factors for the W band are in parenthesis because they are very uncertain (see text). 



20 GHz (AT20G) Survey (|Ekersll2006l ). once completed, will 
cover the whole Southern hemisphere down to 50 mjy. 

Moreover, the Wilkinson Microwave Anisotropy Probe 
(WMAP) mission has produced the first all-sky surveys 
of extragalactic sources at 23, 33. 41, 61 and 94 GHz 
l|Bennett et alj|2003l ; iHinshaw etHlboOTl ). ; From the anal- 
ysis of the first three years survey data the WMAP team 
have o btained a catalogue of 323 extragalactic point sources 
(EPS; Hinsha w et alj (|200?t 0. substantially enlarging the 
first-year one that include d 208 EPS detecte d above a 
flux limit of ~ 0.8-1 Jy |Bennett et alj 120031 '). with an 
estimated completeness limit of ~ 1.2 Jy. The average 
estimated spectral index for detected sources is ~ 0.0, 
thus confirming that the source class which dominates the 
counts at bright fluxes is essentially made of flat-spectrum 
QSOs and BL Lacs, with minor cont ributions coming from 
steep- or inverted-sp e ctrum sources (IToffolatti et al. 1 19981; 
De Zotti et"afl 1 19991 ; IToffolatti et al.l Il999l ; iBennett etaU 



20031 ). 



In order to take full advantage of all the in- 
formation content of the WMAP 3-year survey data, 
lLopez-Caniego et al.l (l200?j. hereafter LC07 ) have used the 
MHW2 filter (|Gonzalez-Nuevo et all 120061 1 to obtain esti- 
mates of (or upper limits on) the flux densities at the WMAP 
frequencies of a complete all-sky sample of 2491 sources at 
\b\ > 5°, brighter than 500 mjy at 5 GHz, or at 1.4 or 
0.84 GHz in regions not covered by 5 GHz su rveys but cov- 
ered by either the NVSS l|Condon et al.lll9981 l or the SUMSS 
(Mauch et al. 2003). This work yielded 5cr detections of 368 
extragalactic sources, including 98 sources not present in 
the WMAP 3-year catalogue. The results were organized in 
a catalogue dubbed NEWPS 5[T (New Extragalactic WMAP 
Point Source) Catalogue. 

As pointed out by LC07, the new flux estimates showed 
small systematic differences, increasing with frequency, with 
those obtained by the WMAP team. These differences were 
attributed to different approximations of the beam shapes. 
In fact, the beams are complex and asymmetrical. Taking 
also into account that sources are observed with different 
beam orientations, it is clear that the effective beam areas 
and, therefore, the flux calibrations, are uncertain. 

To check the calibration of LC07 fluxes we have looked 
for ground based measurements of NEWPSso- sources at fre- 
quencies close to the WMAP ones, finding small but appre- 
ciable systematic differences. This prompted us to investi- 
gate in more detail the calibration problem. The adopted 
methods and the results are presented in §[5] In §[3] we 
present the optical identifications and the redshift distri- 
bution of our sources, based on literature data. After having 



recalibrated the NEWPS5 CT flux densities we have investi- 
gated the spectral index distributions in different frequency 
intervals (§[3J and used them to obtain model independent 
estimates of the mm-wave source counts, beyond the flux 
density intervals where they are directly determined (§0. 
Finally, in §[6l we summarize our main conclusions. 

Throughout this paper we have used the Bayes Cor- 
rected version of the NEWPS Catalo gue to take into ac- 
count the well-known Eddington bias Pddingtonlll94tj| . see 
§ 12.11 or LC07 for more details j]]. The spectral index, a, is 
defined as a e — log (Ss/Si)/log(i/ 2/vi), where Si and 52 
are the fluxes at the frequencies v\ and ^2, respectively (i.e., 
S u oc v~ a ). 



2 STATISTICAL FLUX CALIBRATION 

We have studied in detail the calibration of the estimated 
NEWPS fluxes, after the application of the Bayesian correc- 
tion (see next sub-section) , by adopting two independent ap- 
proaches: a) comparison with ground-based measurements 
at the closest frequencies; b) direct estimates of the effective 
beam areas for the brightest sources. 



2.1 



The Eddington bias and the Bayesian 
correction 



It has long been realized l|Eddingtonl |l94Ch that fluxes of 
sources detected above a certain S/N threshold in noisy 
fields are, on average, overestimated since there are simi- 
lar numbers of positive and negative fluctuations but faint 
sources are more numerous than the bright ones. Thus the 
number of faint sources exceeding the threshold because 
they are on top of positive fluctuations is larger than that 
of brighter sources falling below the threshold because they 
happen to be on top of negative fluctuations. This effect is 
known as the Eddington bias. 

The normalized distribution of true fluxes, S, of extra- 
galactic sources, is usually well described by a power law 

"(1 + 9) 



P(S\q) = kS~ 



S 5? Sn 



(i) 



If the slope q is known, the maximum likelihood estima- 
tor of the true fluxes of the sources is easily calculated 
( Hogg fc Turner! 1 19981 ) from the observed ones. Unfortu- 
nately, in many cases and in particular in the case of the 



1 The Bayes-corrected version of the NEWPS Catalogue can be 
found at |http: / / max.ifca.unican.es/ca niego/NEWPS/ 
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WM AP survey, thi s cond ition is not satisfied. For this rea- 
son, iHerranz" et all (|2006l ) and LC07 have gone beyond the 
results of the previous work by solving simultaneously for 
q and S. We remind here only the asymptotic limits of the 
estimators, valid in the high signal-to-noise regime, refer- 
ring the reader to LC07 for a detailed description of the 
calculations: 



JV 



T7 > In 
N ^ 



(ST) 



Si 



Sf 



1 + 9 



(2) 
(3) 



where is the minimum observed flux and n = S° /o~i is 
the signal to noise ratio of the source. Firstly, we estimate 
the slope, q, from our observed fluxes [eq. @] and then we 
use it to calculate the correction for each source [eq. (3)]. 



2.2 Ground based measurements 

Many WMAP sources are calibrators for the Australia Tele- 
scope Compact Array (ATCA) and/or for the VLA and, 
as such, are continuously monitored at several frequencies. 
A list of 599 sources brighter than 500 mjy at 12 mm (25 
GHz), some of which have also data at 7 and 3 mm (43 
and 100 GHz), is available on the ATCA websitfl We have 
cross-correlated this list with the complete NEWPSso- sub- 
sample flux with S23GHZ 2 Jy. Within a search radius of 
0.348° (the beam-size of the WMAP K-band) we have found 
49 associations, 47 of which have flux measurements also at 
7 mm and 49 at 3 mm. The numbers of NEWPS 5a de- 
tections for these 49 sources are: 48, 48, 43, and 8 at 33, 
41, 61, and 94 GHz, respectively. The median flux density 
ratio for pairs of WMA P channels was obtained using the 
iKaplan fc Meierl (1 19581 ) estimator, taking into account also 
the 5a upper limits. The comparison with calibrators was 
performed in the space of WMAP and calibrator frequencies, 
defining a 5 x 3 matrix of flux density ratios R(vw, Vc), with 
u w = [23, 33, 41, 61, 94] GHz and v c = [25, 43, 100] GHz. 

Since the number of 94 GHz detections is too small, the 
application of the Kaplan-Meier estimator to data at this 
frequency yielded unreliable results. We therefore decided 
not to take into account the _R(94 GHz, u c ) values and to re- 
place them with extrapolations from the other elements of 
the matrix. In practice, for each of the calibrators frequen- 
cies we have linearly fitted the i?([23, 33, 41, 61] GHz,^ c (i)) 
data (where i=l,2,3) and estimated the _R(94 GHz, Vc{i)) el- 
ements by a direct extrapolation of the fits. Next, we es- 
timated the 5x5 matrix R(vw ,vw) by interpolating (we 
have used a 2-dimensional spline interpolation) the esti- 
mated K(i>w,v c ) matrix. The calibration factors for each 
WMAP channel were finally obtained as the diagonal ele- 
ments of the R(yw, vw) matrix. The results and their errors 
are given in the last column of Table [1] Due to the uncer- 
tainty of the correction factor for the W band, hereinafter 
we will use only the corrected fluxes relative to the K, Ka, 
Q and V WMAP bands for our statistical analysis. 



2.3 Self-calibration: methodology 

As stated in § 1, the uncertainty on flux calibration stems 
from our imperfect knowledge of the asymmetric beam re- 
sponse functions. WMAP calibration is based on the bright- 
ness temperature of the CMB dipole; in order to go from 
temperatures to fluxes it is necessary to multiply the tem- 
perature maps by a numerical factor that depends on the 
frequency of observation (thus transforming the maps into 
units of Jy/sr, for example) and then to multiply again 
by the effective area (solid angle) subtended by the point 
source. Therefore, a good knowledge of the effective area of 
the beams is required. In LC07 the effective areas of the dif- 
ferent WMAP beams were calculated using the symmetrized 
beam profiles given by the WMAP tearqj. Obviously, this 
is only a first-order approximation. In particular, for the 
W and V bands the beams are known to be highly non- 
Gaussian and non-symmetric. Thus, accurate flux determi- 
nations require the knowledge of the real, non-symmetric 
effective beam areas at each source position. The WMAP 
team provided real beam maps, but it must be noted that 
the effective beam area at a given position in the sky is the 
result of a complicated composition of the nominal beam 
over many pointings with different orientations at different 
times. 

Fortunately, there is a simpler way to calculate the ef- 
fective beam area, only using WMAP data. In the absence 
of noise of any kind (i.e. if only the point source signal is 
present), it would be sufficient to use known point sources 
to measure directly the beam area: the effective beam area 
would simply be the ratio of the total observed source flux 
to its central value (i.e., the value corresponding to the cen- 
tral pixel). The total NEWPS fluxes have been calculated 
using the formula: 



Ft = A b li(0), 



(4) 



where At is the effective beam area and li(0) is the peak 
intensity for the i-th source, directly calculated from the 

filtered image. 

The SExtractor package (|Bertin &: Arnou"tslll996l ) is a 
well-known software specifically designed for obtaining accu- 
rate photometry of compact astronomical objects on images. 
This software allows us to estimate directly the total flux of 
each source, Fi, without any assumption about the beam 
shape or size. 

To minimize the errors on the estimated fluxes we need 
to select only bright sources, i.e. sources with a high signal- 
to-noise ratio. Therefore, for each WMAP band we have se- 
lected the NEWPS 5(T sources with 1 Jy and local S/N^ 5 
(in the real space) at that frequency; their numbers are 123, 
75, 67, 28, and 3 at 23, 33, 41, 61 and 94 GHz, respectively. 
Then, for each NEWPS source satisfying this selection crite- 
rion we have projected a flat patch around its position. Sub- 
sequently, we have run SExtractor on the resulting patches 
to obtain a direct estimation of the total flux and of its pho- 
tometric error. Besides, for each detected source the peak 
intensity and the background r.m.s. at its position has been 
obtained. Then, we have calculated the effective area of the 



2 http:/ /www. narrabri.atnf.csiro.au/calibrators/ 



3 http:/ /lambda. gsfc.nasa.gov/product/map/dr2/beam_profilcs_ 
get.cfm 
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Figure 1. Redshift distributions of the full \b\ > 10° sample, of 
QSOs and of galaxies (solid, dashed, and dot-dashed histograms, 
respectively). The dotted, dashed, and dot— dashed curves display, 
for co mparison, the predictions of the model by |Pe Zotti et al,l 
||2005|1 for Flat-Spectrum Radio QSOs, BL Lacs, Fanaroff- Riley 
type II sources, while the solid line shows the total. 
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Figure 2. Spectral index distributions for two different frequency 
ranges: 5-23 GHz (upper panel) and 23-61 GHz (lower panel). 
For both distributions the mean, median, standard deviation and 
fraction of flat spectrum sources (see § 14.21) are shown. The ver- 
tical red dashed lines correspond to the median values of the 
distributions. See Section 4 for more details on the calculations. 



beam at the position of the source as Af, i = Fi/U(0). Fi- 
nally, we have computed the average of all the individual ef- 
fective areas, weighted by the corresponding errors, for each 
WMAP channel. 

The results are presented in Table [1] In the first column 
we have the WMAP frequency band. The second column 
shows the average effective beam areas calculated with the 
SExtractor method, and their errors, while the third column 
contains the corresponding beam area calculated using the 
symmetrized beam profiles (this is the beam area used in 
NEWPS for unit conversion) . If the estimated average effec- 
tive areas apply to the whole sky, the correction factor that 
should be applied to NEWPS fluxes is A^sEx/A^sym. This 
factor is given in the fourth column of Table Q] A compar- 
ison with the correction factors estimated from calibrators 
(last column) shows a very encouraging agreement. Since for 
the W band we have only 3 sources that satisfy the selec- 
tion criterion, the correction factor is much more uncertain 
than suggested by the formal error, that only includes the 
dispersion of the estimated values. 

For the following analysis we will correct the NEWPS 
fluxes in LC07 multiplying them by the self-calibration fac- 
tors listed in Table 1, except for the W band, that will not 
be considered further because of the lack of a reliable cal- 
ibration factor, the relatively high number of upper limits 
and its uncertain completeness level. 



3 OPTICAL IDENTIFICATIONS AND 
REDSHIFT DISTRIBUTION 

We have searched the NASA Extragalactic Database 
(NEE0) for optical identifications of the 252 sources at 
|6| > 5° with ^ 5cr detections at 23 GHz and flux den- 
sity above the completeness £?23GHz = 1.1 Jy. The majority 



of sources (181, i.e. 72%) are classified as QSOs, 44 (17%) 
as galaxies, 5 are Galactic objects (2 SNR, 2 HII regions, 
1 PN), and the remaining ones (22, i.e. 9%) are unclassi- 
fied. The galactic objects were obviously dropped from the 
sample. 

The NED database also provided the redshifts of 215 
sources (87% of the extragalactic sample; 95% (42) of galax- 
ies and 96% (173) of QSOs). Although the redshift complete- 
ness is already rather high, it can be increased restricting 
the sample to |6| > 10°. In fact, in the Galactic latitude 
range 5° < |6| < 10° there are 22 extragalactic objects, only 
8 of which have redshift measurements, not surprisingly be- 
cause of the substantial Galactic extinction. Our extragalac- 
tic sample comprises 225 objects at |b| > 10° and 207 (92%) 
of them have redshifts (167 QSOs and 40 galaxies). Opti- 
cal (bj ) magnitudes of these sources were recovered from 
the SuperCOSMOS Sky Survey archive^. The median bj is 
16.7 for galaxies and 18.2 for QSOs. 

The redshift distribution for the latter subsample is 
shown in Fig. [1] where the dashed, dot-dashed, and solid 
histograms refer to QSOs, galaxies, and to the total, re- 
spectively. The median redshift of the sample is 0.86 (0.052 
for Galaxies only, 0.994 for QSOs). We have also plotted, 
for comparison, the redshift distributions predicted by the 
iDe Zotti et all ( 2005) model for different source populations. 
The agreement is generally good, except for the dip in the 
data around log(z) = —1, where the model predicts a little 
bump due to FRII sources. The reason of this discrepancy 
is unclear. The most obvious option is that FRII sources 
are not correctly modeled, but other possibilities, such as a 
large scale inhomogeneity in the distribution of bright radio 
sources, cannot be ruled out. 



4 http:/ /nedwww. ipac.caltech.edu/ 



5 http://www-wfau.roe.ac.uk/sss/ 
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4 SPECTRAL INDEX DISTRIBUTIONS AND 
SPECTRAL CLASSES 

To study the spectral properties of NEWPS sources we de- 
fined a complete sub-sample comprising all the f 9f NEWPS 
sources with 5*23 ^ 1.3 Jy and S/N ^ 5, with counterparts 
in the 5 GHz surveys (GB6 or PMN). We have confined our- 
selves to fluxes somewhat above the 23 GHz completeness 
limit (~ 1.1 Jy) in order to have more than 50% 5a detec- 
tions at 61 GHz (we have a detection rate of 52%); otherwise 
the uncertainties on the spectral index distribution obtained 
from Survival Analysis become unacceptably large. 

In addition to fluxes in the WMAP channels, we have 
included in our analysis 5 GHz measurements. Since the 
angular resolution of the 5 GHz catalogues is much higher 
than that of any of the WMAP channels, we have degraded 
the 5 GHz catalogue to the WMAP resolution. In practice, 
whenever more than one 5 GHz source are present within a 
WMAP resolution element, we have summed up their 5 GHz 
fluxes, weighted with the WMAP beam centered at the po- 
sition of the NEWPS source. To correct for the contribution 
of the background of faint 5 GHz sources we have selected a 
set of 4141 fields (control fields) at 23 GHz with area equal 
to that of a WMAP resolution element at the considered 
frequency and devoid of NEWPS sources. The mean flux 
(120 mjy) of 5 GHz sources found in control fields, weighted 
with the WMAP beam with axis towards the field center, 
has been subtracted from the sum of 5 GHz fluxes (or the 
the flux of the only source) in the fields of NEWPS sources. 
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Figure 3. Examples of typical radiosource spectra at mm wave- 
lengths: a flat-spectrum source (top left panel); a steep-spectrum 
source (bottom left panel); a source whose spectrum flattens at 
v ~ 10 GHz (top right panel); a High Frequency Peaker (HFP) 
source (bottom right panel). Data from the NEWPS Catalogue 
(asterisks) and from the AT20G Survey (courtesy of the AT20G 
group, Massardi et al. (in preparation)) (diamonds) . 



Class 


5-23 GHz 


23-61 GHz 


Flat (%) 


76±6 


59±6 


Steep (%) 


17±1 


41±4 


Inverted (%) 


7±1 






4.1 Spectral index distributions 

The spectral index distributions for the frequency intervals 
5-23 GHz and 23-61 GHz are shown in Fig. [2] The latter 
distribution has been obtained by the Kaplan-Meier (1958) 
estimator, taking into account also the upper limits on the 
fluxes, yielding lower limits on spectral indices. Note that, af- 
ter application of the Kaplan-Meier estimator it is no longer 
possible to distinguish between data and limits. Interest- 
ingly, the distribution of a^-j shows a hint of a second peak 
in the range 0.4-0.8, not seen in the distribution of a| 3 , con- 
sistent with the transition of a subset of sources from flat 
to steep spectra. We have applied the non-parametric two 
sample tests implemented in the public ASURV cod«0 to 
see whether the two distributions are consistent with being 
drawn from the same parent distribution. This hypothesis is 
rejected with a high significance level (more than ~ 8a). As 
noted below, however, the difference may be due, to a large 
extent, to a selection effect. 



4.2 Source spectral classes 

Radio sources show a variety of spectral behaviours in this 
frequency range (cf. I Sadler et al.l 120061) . Figure [3] shows 
some examples: a flat-spectrum (top left panel), a steep- 
spectrum (bottom left panel) , a spectrum flattening at ~ 10 
GHz (top right panel) and a spectrum peaking at high ra- 
dio frequency (bottom right panel). The data are a com- 



http: / /www. astrostatistics.psu.edu/statcodes/sc_censor. html 



Table 2. Distribution of the 191 sources in our sample in three 
spectral classes for the 5-23 GHz and 23-61 GHz frequency inter- 
vals: flat— spectrum (—0.5 a 0.5), steep— spectrum (a > 0.5) 
and inverted-spectrum (a < —0.5). See S I4,2l for more details. 



bination of NEWPS fluxes with the AT20G measurements 
l|Massardi et al.ll2007l '). 

We have classified our sources into three main spec- 
tral classes: flat-spectrum sources (—0.5 ^ a ^ 0.5), steep- 
spectrum sources (a > 0.5) and inverted-spectrum sources 
(a < —0.5). Between 5 and 23 GHz flat-spectrum sources 
dominate the sample: they amount to 76±7% (see Ta- 
bic [H , in agreemen t with conc l usion s of previous work s 
l|Bermett et all 120031 : iRicci et all 120041 : ISadler et all 120061 ). 
However, at higher frequencies, we find a trend towards 
a steepening of the spectral indices. The 5-23 GHz spec- 
tral index distribution has a median (mean) value a = 
0.05 ± 0.04 (0.01 ± 0.03), with a dispersion a = .43, in full 
agreement with the results by IRicci et alj |2004) who used 
the multifrequency ATCA survey data. The median (mean) 
spectral index between 23 and 61 GHz for the same sample, 
selected at 23 GHz, is a = 0.25 ± 0.05 (0.37 ± 0.04) (see 
Fig. [2j ; correspondingly, the fraction of the steep-spectrum 
sources increases from 17 ± 1% to 41 ± 4%. 

The quoted errors are Poissonian in the case of a^ 3 . 
As for of 3, we must take into account that its distribution 
includes a substantial fraction of lower limits, redistributed 
with the Kaplan-Meier (1958) method. Clearly, limits can- 
not be given the same weight as actual measurements in the 
error estimate. As discussed by Cantor (2001), in this case 
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Freq 


Slim 


A (a A ) 


7 7 ) 


Prob 


(GHz) 


(Jy) 


(Jy-isr- 1 ) 






23 


i.i 


34.4 (3.5) 


2.36 (0.10) 


0.84 


33 


l.i 


27.8 (3.2) 


2.25 (0.11) 


0.08 




1.4 


34.6 (5.0) 


2.40 (0.13) 


0.30 


41 


1.1 


30.9 (3.4) 


2.41 (0.11) 


0.99 


61 


1.1 


18.5 (2.8) 


2.24 (0.14) 


0.07 




1.1 


25.8 (4.5) 


2.47 (0.15) 


0.94 




10 10 
log S [Jy] 



10 10 
log S [Jy] 



Table 3. Parameters of the power-law approximation to the dif- 
ferential source number counts calculated from our sample. The 
columns contain: 1. the frequency, 2. the lower flux limits, 3. the 
normalization, A, and 4. the slope, 7, of the counts with their 
errors, and 5. the probability of exceeding the corresponding x 2 
value if the model is correct. Only the data points on source 
counts plotted in Figure 5 and above the lower flux limit, Su m , 
of column 2 are used in this calculation. 




10 10 

log S [Jy] 



10 10 
log S [Jy] 



a simple and conservative estimate of the standard error on 
the redistribution of limits is provided by the Peto formula 



61 
0"23 



N u 



1/2 



(•») 



where P is the probability that a% is larger than the mean 
(or median) value and N u = 91 is the number of lower limits. 
The probability P is approximated by the fraction of sources 
with «23 larger than the mean or the median (P = 0.44 or 
P = 0.50, respectively). The global error is computed as the 
quadratic sum of the Poisson error for the total number of 
data points (measurements plus limits) plus the contribution 
given by eq. ([5]). In the same way, we have estimated the 
errors on the fraction of sources of steep- and flat-spectrum 
sources by identifying the probability P with the fraction 
of sources with spectral index above the boundary for their 
spectral class. 

The 23 GHz selection obviously biases the spectral in- 
dex distribution in favour of flat-/inverted values, which 
yield brighter fluxes at higher frequencies. In the case of a 
Gaussian distribution, a dispersion a a tr anslates into an ap- 
parent flattening of spectral indices by (|Danese fc de Zottil 
1 1984 ICondonll 19841 ): 



Aa ~ (1-7)0-^ ln(23/5), 



(6) 



7 being the slope of differenti al counts. For an Euclidean 
slope (7 = 2.5) and a a = 0.3 (|Ricci et alj|2004 note that 
the distribution in Fig. [5] is broadened by errors on flux 
measurements), Aa ~ —0.2. Thus the observed steepening 
above 23 GHz may be largely a selection effect. 

On the other hand, a real steepening is expected in this 
frequency range since the absorption optical depth decreases 
with increasing frequency at least as v~ 2 so that an increas- 
ing fraction of the emitting regions becomes optically thin, 
i.e. take on spectral indices a ~ 0.7-0.8 or steeper, because 
of electron ageing effects. 



Figure 4. Differential number counts, in bins of Alog(S) = 0.2, 
multiplied by S 2 , of the NEWPSsa sources (circles with Pois- 
son error bars). The so lid curve shows , for co mparison, the total 
counts predicted by the lDe Zotti et al 1 1|2005|) model. The dashed 
lines correspond to the counts extrapolated from 5 GHz, for the 
23 GHz case, and from 23 GHz for the other frequencies, using the 
spectral index distributions of § 14.11 The extrapolated counts are 
too low at 61 GHz, suggesting that the true steepening of spec- 
tral indices above 23 GHz is less pronounced than implied by the 
distribution yielded by the Kaplan-Meier estimator. A good fit is 
obtained damping down the secondary peak of the distribution 
of a%\ (Fig. [2j by a factor of 0.2 (dot-dashed line). 



bv lDe Zotti et alj l|2005l) . The agreement is generally good 
above the completeness limits of the sample. However the 
flux density range probed with sufficient statistics is quite 
limited. 

Over this range, the differential counts can be described 
by the usual power-law approximation: 



dn 
dS 



= AS~ 



(7) 



The best-fit values of the normalization, A, and of the slope, 
7, are given, with their errors, in Table|3j where the goodness 
of the fit is quantified by the probability of exceeding the cor- 
responding x 2 value assuming that the model is correct. The 
normalization, A, is almost constant up to 41 GHz, in keep- 
ing with the fact that the sample is dominated by sources 
with mean spectral index close to 0, but slightly decreases 
at 61 GHz, as the effect of the high-frequency steepening. 

Above 30 GHz, t he values of A are sign ificantly higher 
than the estimates by IWaldram et all (|2007f ). based on ex- 
trapolations of the 9C counts at 15 GHz. As discussed below, 
the discrepancy arises because the spectral index distribu- 
tion of the relatively faint sources of their sample is not 
appropriate at the high flux densities of interest here. 



5 NUMBER COUNTS 

In Fig. [4] we have compared the differential source counts 
of NEWPS 5ct sources at 23, 33, 41, and 61 GHz, normal- 
ized to the Euclidean slope and in bins of Alog(S') = 0.2, 
with the predictions of the cosmological evolution model 



5.1 Number counts estimates by extrapolation 
from lower frequencies 

The spectral index distributions of § 14.11 combined with the 
observational determinations of 5 GHz counts, extending 
over a very broad flux density interval, can be exploited 
to get model independent estimates of the high frequency 
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counts. Alternatively, the comparison of the extrapolated 
counts with the observed ones provides a test of the Kaplan- 
Meier estimates of the spectral index distributions. 

Under the assumption that the spectral index distribu- 
tion is independent of flux density (only valid for limited flux 
density intervals, see below), the source number counts can 
be extrapolated from the frequency v\ to the frequency v 2 , 
convolving the counts at frequency v\ with the probability 
that a source with S\(y\) has Si(y?): 



dN 2 (S 2 (u 2 )) 
d log 62(^2) 



d log 5*2 {v 2 ) 



diVi(Si(i/i)) 
dlog Si(yi) 



P[Si(i/i)|S 2 M]dlog;Si(i/i) 



(8) 



The probability distribution P(Si(vi)\S 2 (v 2 )) is straightfor- 
wardly obtained from the corresponding spectral index dis- 
tribution. 

However, to extrapolate the 5 GHz counts we need such 
distribution for sources selected at the same frequency, while 
that of § 14.11 refers to the 23 GHz selection. In practice, to 
derive the 23 GHz counts knowing those at 5 GHz and our 
distribution of q§ 3 we need to invert eq. ((8)l with v\= 23 
GHz and v 2 = 5 GHz. This is a classical inversion problem 
(Fredholm equation of t he first kind ) and we have used the 
Lucy's iterative method (|Lucvll 19741 ) to solve it: 



dN(S 23 (j)) 
dN(S 23 (j)Y 



r+l 



2- dNr(S 5 (i)) P{S23ij)lS5il)) > 



(9) 



where for the 5 GHz counts, dN°(Ss(i)), we have used the 
description provided by the iDe Zotti et all (|2005r i model, 
that accurately reproduces the data. We have initialized the 
method (i.e for r=0) adopting Euclidean counts, dN/dS oc 
S~ 5 ^ 2 , but we have also checked that our results are insen- 
sitive to the function used to initialize the process. Addi- 
tionally, we have improved the stability of the method by 
smoothing the dN(S 2 3(j)) r+1 in each iteration. 

The results are shown by the dashed line in the up- 
per left-hand panel of Fig. [4] the extrapolated counts are 
in good agreement w ith the direct estimate and with the 
IDe Zotti et al.l (|2005T ) model for S > 1 Jy, but are well above 
the model prediction at fainter fluxes. Since the downturn 
of the co unts below ~ Uy is born out by the 9C counts at 
15 GHz (|Waldram et al.ll2003ft and by th e ATCA pilot sur- 
vey counts at 18 GHz (|Ricci et al.ll2004h . the discrepancy 
must be attributed to a failure of the basic assumption un- 
derlying the extrapolation. This is confirmed by the direct 
evidence of a change with decreasing flux density of the pro- 
portion of flat- and steep-spectrum sources contributing to 
the corresponding 5 GHz counts. This illustrates the risks of 
extrapolations of counts with incomplete information about 
the statistical properties of source populations. 

Additional evidence of substantial changes of the spec- 
tral index distributi on with decreas i ng flu x density at ~ 20 
GHz is provided bv lWaldram et all l|2007l ). who carried out 
mult i- frequency follow-up of a sample complete to 25 mjy 
at 15 GHz of extragalactic sources from the 9C survey. The 
median spectral index between 15 and 43 GHz was found to 
be 0.89 for these sources, much fainter than ours. They also 
reported a spectral steepening with increasing frequency. 

The observed distribution of 023 can be directly ex- 
ploited to extrapolate the 23 GHz counts to higher frequen- 



cies. Since the IDe Zotti et ail (|2005l ) model provides a good 
description of the counts (as determined by the 9C and 
ATCA surveys) also below the limit of WMAP ones at this 
frequency, we have used it to carry out the extrapolations. 
The dashed lines in Fig.|4]show that the extrapolated counts 
tend to be lower than those directly observed, particularly 
at 61 GHz. This is unexpected, since the observed spectral 
index distributions are, if anything, broadened by the con- 
tributions of measurements errors on fluxes and this leads 
to a flattening of the effective spectral index and, corre- 
spondingly, to an overestimate of the high frequency counts, 
contrary to what is found here. 

The likely explanation is that the Kaplan-Meier esti- 
mator over-populates the steep portion of the distribution 
of 0:23. The observed counts are recovered if we scale down 
by a factor of 0.2 the distribution at a 2 ^ > 0.5 (dot-dashed 
line in Fig. [4} . 



6 CONCLUSIONS 

We have presented and discussed a variety of statistical 
properties that characterize the extragalactic sources in the 
NEWPS Catalogue (LC07). The flux calibration has been 
investigated with 2 different methods and correction factors 
have been derived, except for the highest frequency WMAP 
channel, for which the available information is insufficient. 

A search of the NED database has yielded optical identi- 
fications for 89% of the 252 sources in the complete NEWPS 
sample with S23GHZ 1.1 Jy; 5 sources turned out to 
be Galactic, and were dropped. At high Galactic latitudes 
(\b\ > 10°) we obtain a sub-sample of 207 sources (92%) 
with spectroscopic redshifts (again from the NED), for a to- 
tal number of extragalactic sources of 225 . The correspond- 
ing redshift distribution is in gen erally good agreement with 
the predictions of the model by IDe Zotti et all l|2005l ). ex- 
cept possibly around log(z) ~ — 1 where the model predicts 
a significant contribution from FR II sources, that is not 
observed. 

The distribution of spectral indices between 23 and 61 
GHz, Q23j was obtained with the Kaplan-Meier estimator, 
taking into account the upper limits on 61 GHz fluxes. It is 
shifted towards steeper values compared with the distribu- 
tion of alp, obtained using the GB6 or PMN 5 GHz fluxes. 
Although the steepening can be just a selection effect, hints 
of a real high-frequency steepening are noted. 

The source number counts obtained from our sample 
have close to Euclidean slope and are in good agreement 
with the prediction o f the cosmological evolution model by 
IDe Zotti et all i|2005h . 

We have shown that the use of our spectral index dis- 
tribution to extrapolate the 5 GHz counts to high frequen- 
cies leads to inconsistent results below S23GHZ — Uy be- 
cause of a substantial change in the mixture of fiat- and 
steep-spectrum sources. The change, clearly visible compar- 
ing our_distnbjrtonof_spectral indices with that obtained 
bv lWaldram et all ((2007) for a much fainter sample at the 
nearby frequency of 15 GHz, accounts for the difference be- 
tween our high frequency counts and extrapolations from 9C 
counts. This is an example of the risks inherent in extrapo- 
lations with incomplete information of source properties. 

Extrapolations to higher frequencies of the 23 GHz 
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counts using the estimated distribution of a%\ tend to under- 
predict the observed counts, especially at 61 GHz, implying 
that the Kaplan-Meyer estimator probably overpopulates 
the steep portion of the spectral index distribution. The ob- 
served counts are recovered if the distribution is decreased 
by a factor of 0.2 for a%\ > 0.5. 
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